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ABSTRACT 

The threshold voltage for action potential generation is a key regulator of neuronal signal transduction, yet the mechanism 
of its dynamic variation is still not well described. In this paper, we propose that threshold phenomena can be classified as 
parameter thresholds and state thresholds. Voltage thresholds which belong to the state threshold are determined by the 
‘general separatrix’ in state space. We demonstrate that the separatrix generally exists in the state space of neuron models. 
The general form of separatrix was assumed as the function of both states and stimuli and the previously assumed threshold 
evolving equation versus time is naturally deduced from the separatrix. In terms of neuron dynamics, the threshold voltage 
variation, which is affected by different stimuli, is determined by crossing the separatrix at different points in state space. We 
suggest that the separatrix-crossing mechanism in state space is the intrinsic dynamic mechanism for threshold voltages and 
post-stimulus threshold phenomena. These proposals are also systematically verified in example models, three of which have 
analytic separatrices and one is the classic Hodgkin-Huxley model. The separatrix-crossing framework provides an overview of 
the neuronal threshold and will facilitate understanding of the nature of threshold variability. 


Introduction 

Neurons encode stimuli into time sequences of stereotypical membrane potential pulses that are known as action potentials 
(APs). The bring of an AP is thought to be determined by whether the membrane potential exceeds a certain threshold value; 
however, the threshold (following the routine in neuroscience, ‘threshold’ in this paper means ‘threshold voltage/membrane 
potential’ unless particular threshold types are specihed) is not a constant value. Thr eshold variation (i.e., dynamic thresholds) 
has been observed in the electrophysiological experiments both in v/vcP® and in vhrc l^ l ^ l "' 4 Dynamic thresholds can shape the 
responses of neuronPland enhance coincidence detectioiBI ^ as well as improve the feature selectivitjJ^ of single neurons, 
while hlter w eak asynchronous activity in neuronal networkd^liii Threshold variability also participates in precise temporal 
codin^^liZIill and influences metabolic energy efflciencjGSI 

Although the dynamic threshold plays an important role in neuronal information processing, its mechanism has not been 
well characterized. A pioneering study of threshold “accommodation” assumed that the threshold voltage follows a simple 
temporal exponential proces J^. Accordi ngly, similar first-order kinetic equations that describe threshold variation are still used 
by contemporary researcher The equation has the following form: 


de 

Te—= 0„-0, 

dt 


( 1 ) 


where 0 denotes the dynamic threshold, 9oo is a constant threshold and Xg is the time constant. In investigating of the quantitative 
laws of AP generation, Hodgkin and Huxley found that the threshold could be increased by Na^ cha nnel i nactivation and 
channel activation and suggested that the threshold might be a function of the membrane potential 'l l ^" ! Indeed, recent 
experiments have shown that the threshold of cortical neurons in vivo positively correlates with the membrane potential and 
negatively correlates with the rate of me mbrane depolarizatioiPS Additionally, more ion channel types may participate in 
determining the threshold variatiorE3^ In principle, the relationship of thr eshold with either the membrane potential or the 
rising rate is determined by the activation and inactivation of all ion channelJ^lEl!. 

Although the biophysical mechanism has not been well studied by neuroscientists, mathematical classifications and 
mechanisms describing the neuronal threshold, especially the separatrix concept and quasi-threshold phenomena, were 
proposed by FitzHugh in 1950p3. However, these mechanisms have not captured the attention of scientists in a long time, until 
recently, when similar concepts have been hypothesized or reintroduced by researcherP^f^. Several experimental works on the 













threshold variation adopted similar methods or obtained results that were deri ved fro m these concepts; threshold was defined as 
the threshold voltage value at the ti me po int that a stimulus is switched ofiHEH and simple dynamic threshold equations 
similar to equation ([T]) were adopted33iH Additionally, Platkiewicz and Brette deduced a simplified thr eshold equation that 
could quantify the contribution of different ion channels and synaptic conductances to spike thresholdii^li^ and a similar 
form of equation Q is also derivable from their equationJiSl However, the general mechanism of dynamic threshold for AP 
generation in a specified neuron has not be en well described, especially when the neuron is subjected to differing stimuli, 
such as a rectangular pulse or a ramp currenl^lZ^lZZl To connect the theory and experiment results, and considering the r ecent 
advances on the threshold phenomena associated with large difference between the time scales of fast and slow variable£3^ 
that explain the quasi-threshold welP^ we introduce the general notion of sepamtrix, which is a boundary separating two 
different mode s of dyna mic behavior in state space. The general separatrix is equivalent to the concept of ‘threshold manifold’ 
in mathematicand the concept of ‘switching manifold’ in control scienceP^. Two different mechanis ms ex ist that 
describe the threshold voltages in tw o different types of neuronal models: one is the separatrix of fixed pointP^HI and the 
other is the quasi-separatrix or canarcSEHEIIlH The^neral separatrix is able to incorpor ate both the ‘real-separatrix’ (rigid 
mathematical separatrix, the manifolds of fixed pointP^l) and the ‘quasi-separatrix’GHESEH The general separatrix can provide 
an adequate definition for threshold voltage and a uniform mechanism to explain numerous threshold phenomena. 

In the current study, we first propose that the threshold phenomena in excitable neurons can be classified as ‘parameter 
threshold’ and ‘state threshold’. We postulate that threshold voltages belong to the state threshold and have a mechanism 
distinguished from the bifurcation theory that explains parameter threshold. We demonstrate that the general separatrix exists in 
the full state space of different types of excitable neurons. Crossing the separatrix, i.e., separatrix-crossing, at different sites 
causes state threshold phenomena and threshold variation. Then, we introduce a general form of a separatrix in a common 
excitable conductance-based neuron model. According to the form of a separatrix, the general threshold evolving equation 
versus time is deduced and transformed into equation ([T]) for the simplest case. 

Using the separatrix-crossing framework, we then analyze the threshold voltages and threshold phenomena in different 
models. The threshold set (i.e., separatrices) is not a fixed value but is a varying threshold point in a Quadratic Integrate-and-Fire 
(QIF) modeP, a set of threshold curves in a two-dimensional (2D) model, a threshold plane in a three-dimensional (3D) 
piecewise linear (PWL) model and a threshold hypersurface in the Hodgkin-Huxley model. Additionally, the threshold sets of 
the aforementioned models dynamically vary with external stimuli. Through the simple but analytic threshold function in the 
QIF model as well as in the 2D and 3D PWL model, we demonstrate that separatrix-crossing in the state space is one origin of 
threshold variability. We also deduce the threshold variation equations in these models and determine whether the equations 
can be simplified to equation Q. State dependent threshold parameters, e.g., the threshold timing and amplitude of the voltage 
clamp, step current and ramp current, are also determined by separatrix in these models. Additionally, the numerical simulation 
of threshold phenomena in the Hodgkin-Huxley model following voltage clamping demonstrates that the separatrix-determined 
threshold is valid in acmal neurons that exhibit complex dynamics. 

Results 

Classification and mechanisms of threshold phenomena 

Neuronal threshold phenomena can be classified into different types. In general, neuronal threshold phenomena can be 
categorized by their application domains, which are as follows: in vitro, in vivo and in mode^. Threshold phenomena in neural 
models have been mathematically divided into discontinuous-, singular-point- and quasi threshold phenomena (DTP, STP and 
QTP in abbreviation) by FitzHuglP^. In the current study, we propose that the threshold can be generally categorized as either 
the parameter threshold or the state threshold. The parameter threshold depends on bifurcation while the state threshold is 
determined by the initial- or boundary value problem. The two common neuronal threshold phenomena, the threshold voltage 
and threshold current (usually the amplitude of direct current, i.e., DC), belong to the state threshold and parameter threshold, 
respectively. 

The classification is suitable for both theoretical use and practical application. “Parameter” and “state” have distinct 
meanings in models: the parameters will not change, and the states are variables. In reality, a piece of functional nerve 
membrane have several intrinsic statistic properties, for example, the membrane potential, the states of different subunits of ion 
channels, the number of different ion channels available and the reversal potentials for different ion channel types. The former 
two quantities are time-varying, so we call them ‘states’ (or ‘variable’ in the model), while the latter two do not change, so these 
properties are treated as parameters. The biophysical properties in reality have exact correspondents in the conductance-based 
model. Taking the Hodgkin-Huxley model as an example: variables V,m,h, and n form a 4-dimensional state space, while the 
maximum conductance of ion channels g, (/ = Na,K,L) and reversal potential £, [i = Na,K,L) are all parameters. 

In electrophysiological experiments, besides of the intrinsic factors, current can be artificially injected by external circuits or 
received from synapses of upstream neurons. The external current is an independent variable and was treated as a state variable 
by FitzHuglP3. However, the amplimde of DC or step current is unique and can be recognized as a parameter because the 
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current is time-invariant as the parameters do. According to the bifurcation theory, the current amplitude is the most common 
bifurcation parameter. The varying current with special hxed forms, such as periodic or pulse current, may involve parameters 
such as amplitude, duration and period. These quantities may be called ‘parameters’ for the current, but the current varies in 
terms of time. Therefore, the threshold phenomena for these quantities is hybrid. 

Different types of threshold phenomena have different mechanisms. The parameter threshold is usually explained by 
bifurcation theory, e.g., threshold amplitude of DC, while the mechanism of the state threshold is related to the boundary 
separating dynamic behaviors-separatrix. In the current study, we demonstrate that the separatrix determines time-varying 
threshold phenomena, including the threshold voltages and state-dependent threshold quantities, e.g., threshold amplitude and 
time of step or ramp current. 

The classification of parameters and state is rigid in principle; however, considering both the signihcance and time scale of 
variables at the threshold, some approximations can be made. A variable with small significance in threshold variation can 
be approximated as a parameter. The quasi-static effect with large significance should be treated carefully; for example, the 
synaptic current with a long decay may greatly change the behavior of a neuron by slowly varying the threshold in a long time 
scale, which can not be explained otherwise^S These approximations are considered in the modeling process. Once the model 
is established, the classihcation of parameters and states is clear. 


A general mechanism of threshold voltage and threshold variation 

We begin with the general equations of an excitable neuron. Since the AP generation threshold depends on the complex 
interaction of different ion channels and the stimuluJ^DEU, it is necessary to distinguish gating variables, membrane potential 
and external current. On the other hand, considering the important role of adaptation of the ion channels in threshold 
variation, we directly write the time derive of io n cha nnel gating variables in the the adaptation form. In contrast to the 
more general differential equations in mathematic£3l^ we construct the general equations of a single compartment excitable 
conductance-based neuron as 


4 

dt 

= f{v,X)+ie{t), 

(2) 

dX 

X^(v)-X 

(3) 

dt 

Tx(v) ’ 


where v is the membrane potential and X is the vector of the gating variables of the ion channels; they are abbreviations of v(f) 
and X(t), respectively, /(.(f) is the external current stimulus and /(v,X) represents the intrinsic ion currents. The equations 
allow us to focus on the dynamic threshold analysis and assist in deriving the hrst-order threshold equation. 

The threshold voltage of an action potential is not clearly defined, and n eurosc ientists use varying criteria. The most 
commo nly u sed threshold voltage determination is an ad hoc value of dV j (i/1^1212^ obtained by visual inspection. All the 
criteriain^l^ including this method are trying to determine a voltage value standing for the initial point of an AP in the already 
obtained AP time series, and produce some threshold variatiorP^I. To explain and predict the spike-threshold and its variation, 
we define the threshold voltage as instantaneous threshold voltag^^, denoted as 9. The instantaneous threshold voltage is the 
voltage value that, when the membrane potential is instantaneously brought above or below it, will induce one AP or more (see 
Fig. 0 - As opposed to the method in Ref.[^ we do not limit the threshold to be a depolarized value, i.e., a hyperpolarized 
threshold due to inhibition is also acceptable. Theoretically, instantaneously shifting the membrane potential should be realized 
by injecting instantaneous currents (/e(f) = q5{t)), and the instantaneous threshold voltage is searched by adjusting the strength 
of instantaneous current. In the application, a brief current injection is a close approximation of instantaneous current. From the 
dynamic system viewpoint, the instantaneous threshold voltages are dehned by the separatrix, i.e., the boundary that separates 
two behavior modes (the membrane potential returning to the resting state after a small rise or a large excursion). 

As illustrated in Fig. the different separatrix types for different neuron classes define the instantaneous threshold voltages 
well. As an excitable dynamic system, the neuron dehned by equations (|^ and ([^ must possess at lea st one stable equilibrium. 
The stable equilibrium of a neuron is usually a resting state or sometimes a small oscillatiorEinU From the viewpoint of 
a dy namic al system, the resting state is a stable hxed point, and a small amplitude oscillation corresponds to a small limit 
cycid^nm When more than one stable equilibrium exists, the attraction boundary of the resting equilibrium naturally forms the 
separatrix (see Fig|^a, c and d). Usually, the stable manifolds of a saddle (the mathematical real-separatrix) play the threshold 
role and are c alled the STP by FitzHugn^. For those with only one stable equilibrium, the separatrix maybe the canard of a fold 
poini^mmU which is called the quasi-separatri}JSI(pig. |^). The threshold canard always exist s and the threshold phenomenon 
is distinct when the time constants of the predominant fast and slow variables largely diffejum. These phe nomena were 
referred to as QTP by FitzHugh, and can be called canard thresho ld phenom enon (CTP) on the recent advance£311IlH The 
Fenichel’s Theorem ensures the general existence of quasi-separatrijsEHinillllH The typical type I and type II neurons belong to 
the type with more than one equilibrium and the type with only one equilibrium, respectively. However, CTP and STP can 
coexist in a type I neuron (Fig. with enlarged view in c, d and e, e shows CTP). It should be noted that the v-nullcline or the 
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Figure 1. Instantaneous threshold is defined by whether an AP will be induced when only varying the membrane 
potential, (a) Threshold voltages are the minimum/maximum values that determine AP initiation when shifting the membrane 
potential only in state space. Shifting state point s, (in the figure caption, i = 0, 1 or 2) in the v-axis will produce the threshold 
9i in the separatrix. State sq is the resting state, and the injection of a step current pulse with the amplitude I^tep = 0.147 can 
move the state from so to si and S 2 . (b), (c) and (d) Time series of membrane potential (v) and external current (ie). 
Instantaneous threshold 0, determines whether the injection of a brief current will cause state s, to generate an action potential. 
The durations of the step current injection to si and S 2 are 7.429 and 12.5, respectively. The corresponding thresholds in (b), (c) 
and (d) are also plotted (red dotted lines) and marked. The model is a classic FHN model with the following parameters: 

= 15, a = 1.25,17 = 0.875. 


curve that dv/dt equals a certain values is not the threshold voltage set, especially in Fig. in which both initial values exhibit 
a negative dv/dt. The threshold voltage determined by an ad hoc large positive value of dv/dt is the indication that an AP has 
happened but is not a predictive criterion of whether an AP will be generated. 

As described in the above paragraph, a separatrix is a set of special state points that determines the threshold for AP 
generation in the phase/state space. For a certain neuron, a separatrix is the function of the state variables and external current, 
but the specific form depends on parameters of channel gating variables X, such as maximum conductance gx and reversal 
potential The full expression of the separatrix should be written as S (v,X; 4; A) =0, where A denotes the vector of 
parameters. Focusing on the dynamic threshold of a specific neuron, we ignore the specific form of a separatrix and assume that 
a separatrix has a general form of S{v,X',ie) = 0, so the voltage values at the separatrix, 9{X',ie), are instantaneous threshold 
voltages. When the state of a neuron is changed by a stimulus from one side to the other of a separatrix in state space, the neuron 
switches between AP generation and subthreshold fluctuation. The crossing at different points of a dynamic separatrix caused 
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Figure 2. Separatrices determine threshold voltages in the state space of type I and type II neurons. Red and magenta 
solid lines are separatrices. The black solid lines are real trajectories, and the dashed lines are nullclines (dv/dt = 0 and 
dw/dt = 0), respectively, (a) Separatrices of type I neuron; Boltzmann-FHN neuron model. The red separatrices represent the 
two stable manifolds of the saddle (left open circle), one extends to minus infinity while the other connects to the unstable node 
(right open circle). The magenta separatrix is in fact a trajectory that connects the unstable node and the fold point of cubic 
v-nullcline. Rectangles indicate the enlarged regions of (c), (d) and (e). (b) Separatrix of type II neuron; classic FHN neuron 
model. There exists only one equilibrium that represents the resting state, and the canard-mechanisrrPH determines the 
threshold set (quasi-separatrix) of this type II neuron, (c), (d) and (e) The enlarged regions of (a) show the coexistence of 
separatrices and the separatrices determine the thresholds rather than the v-nullcline or the curve in which dv/dt is equal to a 
certain value. Up/left and down/right blue dashed lines represent the curve that dv/dt equals —0.015 and 0.015, respectively. 


by different stimuli bring about dynamic spike thresholds, so the separatrices of a neuron determine the possible variation scale 
of the spike threshold. 

According to the general form of a separatrix, i.e., 0(X;ie), and the normal form of the time evolving equation of ion 
channel gating variables (eq.([^), we observe that the threshold varying versus time; 


dQ (X, ie) 
dt 


do die{t) 

die dt 


a:GX 


do dx 
dx dt 


do die{t) 

die dt 


+ E 

j:GX 


do Xoc(y)—x 
dx X^iy) 


(4) 


This equation indicates that the dynamic threshold depends on both the external current and the adaptation of all of ion 
channels, in principle. For neurons with DC injection or without external inputs, the first term in the right hand of the above 
equation is zero. Additionally, if one channel gating variable dominates the threshold dynamics, then equation Q can be 
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simplified as: 


dd dd 


(5) 


If do jdx is a constant, i.e., thresholds are linearly dependent on the dominant gating variable (9(X) ^ kx, here A: is a constant), 
and the characteristic time (t^(v)) is approximately a constant for sub-threshold currents, the above equation can be further 
transformed into the simplest threshold dynamic equation-equationQ. Equation Q is directly assumed in Refs.[9][T5|an 
while in Ref. 16 it is deduced by a special condition that has no simple correspondence. Our derivation implies that equation 
(1) will fit well for the situations in which the threshold is linearly varied with the only one dominant modulating variable with 
a approximately constant characteristic time, in the process of DC injection or when no external current is presented. 

Threshold states explained by separatrices are distinctly distinguished from and, at the same time, are associated with 
threshold parameters described by the bifurcation theory. The bifurcation theory describes different dynamic behavior caused 
by the variation of fixed parameters and does not take into account the initial states. The separatrix determines the threshold 
phenomenon of different initial values. Take the external current injection as an example; the bifurcation theory describes the 
threshold amplitude of DC and relate to the rate coding of periodic repeated firing, while the separatrix framework explains the 
threshold voltages under subthreshold DC or transient inputs (e.g., step current) and explains the temporal coding. However, 
parameters determine the specific form of a separatrix. The variation of parameters will quantitatively change the separatrix if 
no bifurcation occurs. Parameter thresholds determined by bifurcation change the form of the s epara trix, and may even cause 
the separatrix to appear or disappear, e.g., saddle-node bifurcation in the type I neuron modeEilH When entering into the 
mode of repeated firing controlled by a global limited cycle, separatrices/state thresholds disappear. To completely predict the 
response of a neuron to the varying stimuli, we should combine both mechanisms. 

Recent experiments and numerical simulations have attempted to elucidate the i nfluence of the external stimulus on 
threshold voltages by defining threshold voltages as the post-stimulus threshold valueJ^HEH The post-stimulus threshold 
is a unique case in that ig{t > fg) = 0 (4 denote the time that the stimulus is off), i.e., the separatrix after the stimulus off is 
0(X;O). The separatrix is therefore the same as the old separatrix prior to the stimulus. The mechanism for post-stimulus 
threshold indicates that the changed state crosses the same separatri x. Sep aratrix-crossing in transient responses is different: the 
almost unaltered state left for the other side of the varied separatri5(E2lII[ We emphasize the word ‘almost’, because the abrupt 
injected current does change the state; however, the magnitude of voltage shifting is far smaller than the shifting of v-nullcline 
(/At/C ^ /), when the switching duration Af is very short. In more general situations, the state and separatrix vary at the mean 
time with the continuous varying stimulus. 

A voltage clamp Axes the membrane potential, but it changes other states of the membrane and therefore also alters the 
threshold. After a voltage clamp, AP generation is also determined by whether the clamped voltage exceeds the instantaneous 
threshold 0(X;O). After any stimulus, if no AP is produced following the most recent deviation from the resting potential, then 
whether the system will produce an AP is determined by whether the state after the stimulus is turned off crosses the separatrix 
0(X;/g). 

In the following paragraphs, we apply the separatrix to different models and explain threshold phenomena using the 
separatrix-crossing mechanism. The separatrices, i.e., threshold sets, of these model are a point in the QIF model, a curve in 2D 
models, a plane or surface in 3D models and a hypersurface in the four dimensional HH model. 


Dynamic threshold point determined by stimulus 

In the classic Integrate-and-Fire (IF) model, the threshold voltage is a fixed parameter . Rece ntly, the dynamic threshold has been 
introduced using artificial mechanisms to predict the spiking series of real neuron^l^E^ Here, we demonstrate the dynamic 
threshold in a one-dimensional quadratic integrate-and-flre (QIF) modePil. The separatrix in a QIF model is a point, which is 
limited by parameters but varies in response to varying stimulus strengths. The dimensionless equation of the QIF model ipa 

= {v-Vr){v-Vt) + ie{,t), if V > then V ^ v^scr- (6) 

The parameters are related in that Vp^a^ > v, > Vr > Vreset- When ig{t) = 0, the resting potential is Vy and the threshold is v,. 

According to the bifurcation theory, the QIF model is topologically equal to the canonical type I neuron model and exhibits 
a saddle-node bifurcation if the DC current amplitude crosses the threshold value Irehobase- hehobase is determined by the 
equation [vt + VrY —A[vrVt+Ie) = 0, i.e., the discriminant of the square equation dv/dt = 0. For a DC input with the amplitude 
4, if 4 > Irehobase, the neuron periodically discharges; if 4 < Irehobase, the two real, unequal roots of equation dv/dt — 0 are 
two equilibria, threshold 0 and resting potential Vrest- 

For the varying external current ie{t), replacing 4 with ie{t) leads to the dynamic threshold 0(ie) and dynamic resting 
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potential Vresti^e) of- 




Vt + Vr± y/jvt + Vr)^ -^jVrVt +ieit)) 
2 


(7) 


Both the threshold 9{ie) and resting potential Vrest{ie) vary with the external current ie{t). The subthreshold depolarizing current 
elevates the resting potential and lowers the threshold, i.e., decrease the difference between the threshold and resting potential. 
The hyperpolarized current increases the difference between 9{ie) and Vrest{ie), vice versa. When parameters are fixed in the 
QIF model, the threshold voltage is only related to the external stimulus. As displayed in Fig. [^, the same membrane potential 
may be subthreshold with no stimulus or suprathreshold when a sufficiently strong current was injected. Once the membrane 
potential crosses the threshold, the neuron is in the process of AP generation. AP generation is a fast, but not instantaneous, 
process, and can be interrupted even during the depolarizing phase. After the moment the threshold is crossed, a strong and 
long enough hyperpolarized current injection, which elevates the threshold, may interrupt the firing process (see Fig. [^). In a 
rigid sense, after the membrane potential exceeds the instantaneous threshold, the application of external current ensured a 
persistent increase in the difference between the membrane potential and the instantaneous thresholds determine the AP. To 
predict the generation of an AP, it is necessary to know both the subsequent stimulus and the present state of the neuron. 
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Figure 3. Stimulus determine the threshold point in the Quadratic Integrate-and-Fire model, (a) After a ramp current, 
an AP cannot be induced without stimulation (blue solid line); a subthreshold constant current that lower the threshold can 
elicit an AP (red solid line). The different behavior can be predicted by the threshold lines-corresponding dashed lines. Red 
lines represent the situation with subthreshold current Ig = 50. The black dash-dotted line represents the threshold current, (b) 
The rectangular pulse induced AP (blue solid line) can be prolonged (red solid line) or prevented (green solid line) by 
hyperpolarized current injection. Whether the membrane potential exceeds the corresponding threshold predicts whether the 
neuron will generate an action potential after the stimulus switching. 


Dynamic threshold curve in two-dimensional models 

As a one-dimensional system, the QIF model has a threshold value that varies with an external stimulus. In the similar 
2-dimensionaI Izhikevich model, the threshold expands to a curv^Sl xhe threshold cur ves in a 2-dimensionaI model have 
been demonstrated in many dynamic syste ms, suc h as the FitzHugh-Nagumo (FHN) model2§12111i|, the reduced 2-dimensional 
HH modeiSI the Morris-Lecar (ML) modeEllll, etc .EH However, the quantitative form of the separatrix cannot be precisely 
determined in these models. We present the analytic separatrix in a bi-dimensional model with piecewise linear nullclines and 
explain, in detail, the threshold phenomena that are determined by the separatrix-crossing mechanism for voltage clamping and 
other stimuli. We also deduce the first-order differential equation that describes the dynamic threshold variation versus time. 
This 2D PWL model is analytically solvable and can be regarded as a modified FHN model. 

As shown in Fig. the real-separatrix determine the thresholds of our 2D PWL model. The real-separatrix is a straight line 
that is precisely determined by the intersecting point (saddle) of the middle segment of thev-nullcline and vr-nullcline. The 
separatrix in the hyperpolarized voltage region is in fact the trajectory connected to the real separatrix. 
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Figure 4. Separatrices and the separatrix-crossing mechanism after voltage clamp in the state plane of a 2D PWL 
model. Blue and cyan dashed straight-lines represent the v-nullcline and w-nullcline, respectively. Black dots and circles are 
three fixed points; the left one (stable spiral, black dot) is real, while the middle one (stable spiral, black dot) and the right one 
(saddle, black circle) are virtual. Wide red dash-dotted curves delineate the global separatrices, the right straight one in the 
middle region is a real separatrix (stable manifold) of the saddle, while the left winding one is a trajectory connected to the real 
separatrix. Blue and green solid lines are real trajectories. Black lines show the voltage clamp process, the dashed one 
demonstrates the ideal instantaneous voltage shift and the solid two with an arrow marking the time evolution show the voltage 
holding after the sudden shift. Initialized at state points 1,3,5 or 7, the neuron will follow the blue trajectory and has a small 
peak, which is considered a subthreshold fluctuation; initialized at the state points 2, 4, 6 or 8 (on the other side of 
separatrix-firing region), the peak of the membrane potential will be large and will be considered an AP. The difference 
between 1 and 2 (5 and 6) is whether the sudden shift of membrane potential cross the separatrix, while the difference between 
3 and 4 (7 and 8) is whether the continuous voltage holding crosses the separatrix. 


The real separatrix of our 2D PWL model is a straight line in the v — w state plane 
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where k„, k„, and bm are parameters; and kQ is the slope rate of the separatrix which depends on the parameters including the 
time constant Xw (see Section Methods). The threshold is linearly related to the recovery variable w. Because kg > km> k„ 
(see Section Methods), we can infer that the threshold decreases as the strength of the depolarizing current increases. The ig 
increment linearly moves the separatrix (a straight line) up in the v — w plane. Therefore, an external current linearly changes 
the thresholds. 

In the following paragraphs of this subsection, we will demonstrate that separatrix-crossing determines threshold parameters 
and threshold voltage variation after voltage clamping, step current and ramp current injections. 

We first demonstrate the separatrix-crossing mechanism of threshold variation after voltage clamp (0(w, 0)) in the state 
plane of our 2D PWL model in Fig. For a clear demonstration, we only use two different trajectories that originate from the 
state points 1 and 2 in the phase plane with the former (blue solid line) on the non-firing zone and the latter (green solid line) on 
the other side of the separatrix (red broad dash-dotted line) representing the bring zone. From Fig. we observe the threshold 
phenomenon in rapid voltage shift from resting potential: the instantaneous voltage shifts from state 1 to state 2 (from state 
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5 to state 6) horizontally across the separatrix from the resting zone to the firing region and thus an AP was generated. The 
continuous voltage holding vertically evolves and then crosses the separatrix in the state plane, thereby determines another 
threshold phenomenon: the threshold timing of a voltage clamp at a fixed voltage. For example, if a depolarized voltage is 
held at ISmV and the holding time is lengthened, the neuron will pass through state point 7 to state point 8, which prevent AP 
generation. Inversely, if a hyperpolarized voltage is held at — l5mV, an increased clamping time will cause an AP generation 
(state 4) while a shorter clamping will not (state 3). The long voltage clamping at hyperpolarized potential facilitates AP 
generation. Therefore, the winding separatrix around the resting state in the hyperpolarized region results in post-inhibitory 
facilitation. 

As shown in Fig. the threshold voltages and timing that are determined by separatrix-crossing are clearly exhibited by 
the maximum depolarizing voltages after voltage clamps with clamped voltage Vc and clamping duration T^. Comparing two 
subhgures in Fig. we find that a long voltage holding at a depolarized potential hampers AP generation. Conversely, a long 
voltage holding at a hyperpolarized potential facilitates AP generation. This phenomenon is caused by separatrix-crossing 
from the firing zone to the subthreshold zone during depolarized voltage holding and vice versa during hyperpolarized voltage 
holding, as demonstrated in Fig. |^(from state 8 to state 7 and from state 4 to state 3, respectively). The border that is indicated 
by the color in the low voltage region (v < 25mV) in Fig. represents the normally mentioned threshold voltage and correctly 
fits the curve Vc = 9(vc, Tc), which is the projection of the real-separatrix in the Tc — Vc phase plane. The border in the high 
voltage region corresponds to the line dv/dt = 0. Additionally, in Fig. |^, the separatrix projection in the hyperpolarized 
voltage region (left region in phase space) precisely limits the domain of the bring region of the numerical results. 
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Figure 5. Maximum voltages indicate threshold clamping voltages and times in the 2D PWL model. Numerical results 
(colored patches) fit the analytical results (colored lines), (a) Maximum voltages after depolarization voltage clamp. For a fixed 
clamping voltage Vc, threshold clamping duration (t^) in the low voltage region (analytical results are indicated by the magenta 
dash-dotted line) indicate that a long depolarized voltage suppresses AP bring, while threshold clamping duration for a high 
voltage well bt the curve dv/dt —Q (yellow lines) which indicates that longer clamping potential time-series is the lack of a 
depolarizing phase of AP. (b) Maximum voltages of hyperpolarized voltage clamping. For a constant voltage, inverse with the 
depolarized situation, a short voltage clamping cannot induce an AP but a longer voltage clamping can. The threshold voltage 
for a bxed clamping time is a monotonic increasing function of clamping duration. 


Similar to voltage clamping, current injection, such as rectangular and ramp current pulses, will push the state across the 
separatrix if the pulse duration is long enough (see Supplementary Fig. SI). The threshold pulse strengths and durations, which 
are determined by real-separatrices in the v — w phase plane, may form a threshold boundary in the parameter plane. 
According to eq.Q and ([^, we can deduce the brst order threshold equation 


dQ 

dt 


9oo{w,ie) -9- 


^6 i]^m 


(9) 


If dig/dt = 0, i.e., 4(f) = 4, our 2D PWL model satisfy the conditions of equation Q, and we obtain the same threshold 
equation reported in Ref. an dED 

Our PWL model is more like a type III neuron according to Flodgkin’s classibcation, as it can not generate periodic APs. 
However, if we set the slope of w-nullcline to be larger than the middle segment of v-nullcline, the PWL model is capable of 
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firing periodically and of having a type II /-/ curve. Additionally, the separatrix is an unstable manifold of the unstable node 
of nullclines in the middle region (See Fig. S2.). Our 2D PWL models have only one real equilibrium and lack fold points 
(v-nullcline where has dw/dv = 0) between the middle segment and right segment of v-nullcline, and more important is the 
distinctiveness of threshold phenomena depends on the time constant (see Fig. S3.). The threshold phenomena in these models 
essentially belong to CTP (QTP). However, it is interesting that the v-nullcline is discontinuous (DTP), and the threshold can be 
expressed as the real separatrix which locates in the layer of canards (STP). 

Threshold surface in three-dimensional neuron models 

Here, we extend the two-dimensional PWL model to a three-dimensional model, which provides an example of three- 
dimensional phase space analysis of threshold phenomena. 
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Figure 6. Separatrix determine the threshold voltages in three-dimensional PWL model. Trajectories 1, 3 and 5 (blue, 
red and magenta lines) behind the threshold plane and trajectories 2, 4 and 6 (green, cyan and yellow lines) before the threshold 
plane exhibit small difference in initial values; however, they have great behavior difference, subthreshold fluctuation or firing. 
The gray transparent plane, which divide the space into subthreshold and suprathreshold regions, is a separatrix determined by 
stable (red arrow) and unstable (blue arrow) eigenvectors of an equilibrium (blue point). 


As shown in Fig. a plane of threshold states that separates the 3-dimensional space into two parts defines the firing 
or non-firing region in the middle region. The separatrix in this 3D PWL model is a real separatrix that is actually a plane 
determined by the stable and unstable eigenvectors. For dynamic analysis, see Supplementary Note SI and Supplementary Fig. 
S4. 

The separatrix, i.e., a threshold plane in the middle region is 


0(m,w, L) = —au — bw + 


(1 +aku+bk„) 

k +k -k + 


( 10 ) 


where a, b, ku, k„ and k^ are parameters (see Methods section). The threshold variation equation versus time can be obtained 
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as 


dO kuV — u ^k„v — w Iaku+bk„ dig 

— =—a - b -h ■■ —. (11) 

dt Xif Xw kii “t" kyii dt 

The two characteristic time constants, T„ and T»., prevent the simplification of equation into equation Q if the difference 
of Th and x„ is neither very large nor very small. 

The threshold variation phenomena of post-stimulus, such as a voltage clamp, step current or ramp current, are also 
determined by separatrix-crossing in the 3D PWL model. Additionally, it is a natural extrapolation that similar higher 
dimensional PWL models may have separatrices in the form of hyperplanes. 

Similar to our 2D PWL model, the threshold is in this 3D PWL model is essentially belong to CTP, but the threshold set can 
be expressed as the real-separatrix which is in the layer of canards. As shown in Figure Sx, we also provide the continuous 
threshold hypersurface in 3D FHN model, which is determined by a quasi-separatrix in 3D state space. The fold point in Fig|^ 
expands to a fold curve (see Fig. S5, the straight line ). The separatrix comprises threshold manifolds passing through the fold 
points in the fold curve. 


Threshold hypersurface in the HH model 

We numerically demonstrate the separatrix-crossing mechanism in the classic HH model after voltage clamps, as shown in Fig. 
1^ The maximum depolarizing voltage V,nax (indicated by color) after a voltage clamp under the corresponding Vg (clamped 
voltage) and Xg (clamping duration) is mapped onto the voltage clamping surface. The clear boundary of the maximum voltage 
is a numerical approximation of the separatrix, which shows distinct characteristics in subthreshold and firing regions. The 
similar separatrices can be numerically plotted in other 3D state spaces (such as the m-h-n space, see Supplementary Fig. S6). 
The separatrix line in the voltage-clamp surface is part of a high-dimensional separatrix that intersects the voltage clamping 
surface. 

Voltage clamp transient processes at several different potentials are displayed as yellow solid curves with the potential 
marked at the end (Vc-lines) in Fig. while the green dashed lines with time marks (Tc-lines) show the voltage clamping 
duration. The red dash-dotted line indicates the asymptotic value as time approaches infinity under different clamp voltages 
and h^{V)). As Vc-lines indicate, when the voltage is clamped at a certain value, an immediate release causes the 
system to remain in the firing zone, allowing for the generation of an AP. Conversely, if being held at a voltage for a long time, 
the system will cross the separatrix into the subthreshold zone, and APs cannot be generated. 

As seen in the PWL models, whether the HH neuron will generate an action potential after voltage clamping is determined 
by whether the clamped voltage (V^) exceed the instantaneous threshold, 9. We numerically show the scenarios of threshold 
variation during several voltage clampings with varying clamped voltages and clamping durations in Supporting Information 
(Supplementary Fig. S7). The intersecting points of Vg and 9 form the clear threshold border in the low voltage region of 
Fig.li.e., Vg = 9{Vg,Xg',0) as m, h and n is the function of V^ and Xg according to equation ([29])). The similarity of maximum 
voltages versus clamping voltages and durations in the 2D PWL model and the HH model indicate similar separatrix-crossing 
mechanism, but the differences in the details indicate the complex interaction of different gating of ion channels in the HH 
model (Supplementary Fig. S8). 

According to the general equations Q and ( |29| ), the instantaneous threshold voltage variation in the HH model following a 
voltage clamp {9{Vg, Te;0)) is determined by 


d9{m,h,n) 

dt 


y 99 1 

“ dx xyVg) 


{xoo(Vg)-xo)e 


(x = m,h,n) 


( 12 ) 


From the above equation, we observe that in principle, all channels and their subunits exert effects on the threshold variations 
because they are all coupled to the membrane potential; however, the significance is determined by weights of both the time 
constant Xx{Vg) and d9 jdx. Additionally, if the external stimulus is not a voltage clamp, e.g., ramp current, the ‘time constant’ 
Xxiy) for each gating variable also varies with the membrane potential, and therefore, the temporal evolution of post-stimulus 
threshold is more complex. 


Discussion 

Highly variable spiking thresholds are observed both in and in vivcP®! Recent studies concerning the threshold 

variation of neural spikes demonstrate three distinct charact eristics. First, as the threshold is stimulus dependent, most 
researchers use the threshold after the stimulus is ofl®31Zl In the current study, we show that the separatrix-crossing 
mechanism is universal for post-stimulus threshold phenome na. Se cond, because the threshold is not a single value, several 
studies have extended the threshold point to a threshold curv£3^ In fact, the separatrix concept has been proposed since 
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Figure 7. The separatrix-crossing mechanism of the HH model after voltage clamp in the V-m-h state space. Colors 
on the voltage clamping surface indicate maximum voltages after voltage clamping. The clear boundary between blue and other 
colors is the numerical approximation of separatrix on a voltage clamping surface in the V-m-h phase space. The red 
dash-dotted line reveals the asymptote when the voltage clamping time is infinity. Yellow solid lines show voltage clamp 
processes at several voltages (Vc-lines) and the green dashed lines indicate the same duration after voltage clamp (Tc-lines); all 
are marked with a number, showing corresponding values at the end. The separatrix in V-m-h space is part of the global 
separatrix in four dimensional phase space and divides the voltage clamping surface into firing and no firing regions A Vc-line 
shows that the state of the neuron continuously changes, crosses the separatrix and comes into the no firing region as the 
clamping time Tc increases. 


the 195 o£ 3EH and similar hypotheses for general threshold also have recently been reintroducecP^^^. In the current paper, 
we demonstrate that different types of separatrices generally exist in excitable neuron models and systematically show that 
separatrix-crossing mechanism determines threshold phenomena in example models, especially the 1-, 2- and 3-dimensional 
models with analytic separatrices. We also found that the singular-point threshold and canard threshold phenomena coexist in 
type I neuron which have previously been considered to only have singular-point thresho ld phenom ena. Third, researchers 
attempt to simplify the threshold problem and construct quantitative threshold equation^E^nilH Incorporating dynamic 
thresholds into the IF model was shown to be a very efficient way to predic t spik e events in real neuronJSl Additionally, 
researchers have simplified the threshold variation to a threshold equatiorU^EEIl Here , we deduc ed a detailed threshold 
evolving equation versus time, which can be simplified to the previously assumed equatiorSHESGlIlIlfor the simplest case. 

The threshold voltages were positively correlated wit h the m embrane potential (vjDS and negatively correlated with the 
rising rate of membrane potential (dV/dt) before the spikJ35Ell. Our present framework is deterministic, we discuss these two 
relationships in identical deterministic processes. The mean potential prior to the spontaneous spikes approximate the resting 
potential. The mean potential measurement is a long enough process (250ms in Refs.[T]-[^, so we can assume that the state of 
the neuron in different mean membrane potentials is the resting state on that membrane potential. Such a process is approached 
by voltage clamping, and the final resting states of different membrane potential approach the w-nullcline. Shown in Fig. S9, 
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Figure 8. Maximum voltages of HH neurons indicate threshold clamping voltages and times, (a) Maximum voltages 
after depolarized voltage clamp. For a fixed clamping voltage 14, as the clamping duration increases, the peak of the action 
potential decreases. Extended voltage holding inhibits AP generation after voltage clamping. The border in the low voltage 
region displays the normal threshold and the border in the high voltage region is determined by dV (dt — 0 (see Supplementary 
Fig. S8). (b) Maximum voltages of hyperpolarized voltage clamps. For a constant voltage, inverse with the depolarized 
situation; the short time clamps cannot induce an AP while the long-enough ones can. As opposed to the depolarized situation, 
the threshold voltage for a fixed clamping time is a monotonic increasing function. 


the thresholds in w-nullclines do increase with membrane potential. The relationship of the threshold to the rising rate of the 
membrane potential is more obvious. As demonstrated in Fig. S9, in the Boltzman-FHN model, a stronger step current with a 
larger depolarizing rate (dv/dt) produces a less depolarized threshold. As displayed in both Fig. SI and S9, the stronger step 
current cause a lower trajectory which means the recovery variable w varied less for the same v increment than that of the weak 
current. The stronger the current, the lower the trajectory. The infinite strong current injections shift the membrane potential 
and hold the recovery variable, which is the ideal current form that defines the instantaneous threshold voltage (from state point 
so to Bo). 


Our framework includes different ion channel gating variables with different characteristi c time s. The ion channel subunits 
with long characteristic times may induce the threshold variation related to previous spikeJ^lEH According to Ref. 


34 


the 

synaptic current with a long decay can determine the cessation of rebound spiking, which can not be explained by inhibitory 
current, as if the synaptic current was added using a differential equation of an independent variable. Thus, we infer that 
the separatrix in full state space should include the threshold modulation mechanisms over a long time, if all of the factors 
influencing the threshold are included in the differential equations. If the difference in characteristic time scales between 
dominant slow and fast variables is not large, intermediate-amplitude AP will be generated and threshold phenomena maybe 
ambig uouJSI w hich is the the reason why the quasi-separatrix is always mentioned as ‘not well-defined threshold’ in neuron 
modelj^l ^^l^^ l However, considering that the type II neurons are functional in nervous systems, so it is reasonable to infer 
that either these neurons do have large enough differences of characteristic time scales, or other factors such as noise are able 
to overcome this weakness. In addition, the threshold voltages can still be defined by the quasi-separatrix even though the 
threshold phenomena are not distinct. 


In the present framework of separatrix-crossing, the noise or highly fluctuating inputs are not included. In our framework, 
the noise and the fluctuating inputs do not only change the state of the neuron but also vary the threshold. However, because 
that the threshold is instantaneous, it is possible to extend the framework to including the noise or fluctuating inputs, which 
cause instantaneous random variation of the threshold and state. Khovanov et. al. have described how the external and internal 
noise cause the spontaneous spikes in the FHN model by considering the quasi-separatrijJ^. 

The threshold variability and the rapidity of onset of AP have inspired a rece nt deb atd ^l^l^^ EHEH concerning whether the 
two characteristics can be satisfied at the same time by using th e clas sic HH modeC®3. Two explanations that try to reconcile 
the two conflicting facts concern more than one compartmeni^H^^ In a multi-compartment model, the thresholds of each 
compartment are also affected by the neighbor’s states, such as the lateral currenP^and the relative size of the soma and axorP^ 
Evidence from experimental, theoretical and numerical simulations demonstrate that the morphology and distribution of ion 
channels can affect the thresholds of a neuron, especially the axon initial segment, which has a very high sodium channel 
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densitjGnmilH How will the separatrices of compartments determine the threshold of the neuron? Our simple PWL models 
with precise thresholds provide applicable models that are suitable for future exploration of these concepts. 

In summary, by introducing the general separatrix concept, we demonstrate that the separatrix generally exists in excitable 
neurons, and determine threshold voltages, threshold voltage variation and post-stimulus threshold parameters for AP generation. 
We brought forward the idea that threshold phenomena can be classified into two general types with different mechanisms, 
the parameter threshold, described by the bifurcation theory, and the state threshold, explained by separatrix-crossing. The 
separatrix-crossing mechanism also determines the the threshold of timing-related parameters of independent variables. We 
also natur ally de duced the threshold variation equation versus time and transformed it into the simple form used by previous 
researcher jmmHIll for a simpler case. Separatrix and separatrix-crossing-determined threshold phenomena were demonstrated 
and verified in example models. The separatrix-crossing mechanism and models we proposed will benefit further investigation 
of threshold variability on other factors, such as morphology, passive membrane properties and the spatial distribution of the 
ion channels of neurons. 

Methods 

Classic FHN model and Boltzmann-FHN model 

The classic FHN model is described by differential equations; 


dv 



dt 

= V- — -W + le{t), 

(13) 

dw 

kwV + b„ — w 

(14) 

dt 

Tvp 


where v represents membrane potential, and w is a recovery variable and ig(f) is time-varying injected external current. To 
demonstrate of canard threshold phenomenon well, we let = 15, k„ = 1.25, b„ = 0.875 and ie{t) = 0. The neuron belongs 
to type II neuron. 

Similar to ML-FHN neuron model used in Ref. we modified the equation of recovery variable w to: 
(l+exp(-f7(v-c)) 

As the nullcline of w is a Boltzmann function, we call it as ‘Boltzmann-FHN model’. To show three separatrices determined by 
saddle and canard phenomenon well, we take T„ = %, a = 2,b — h, c = 0.27 and 4(f) = 0.62. The model is a type I neuron. 

To show how canard trajectories constitute a quasi-separatrix in three dimensional state space (Fig. S5), we also extend the 
classic FHN model to be three dimensional by adding a similar slow recovery variable u with the same form of equation ( [T4l l. 
The parameters for 3D FHN model are as follow: = 25, k„ = 1.25 and b„ = 0; = 15, ku = 1.25 and = 0; / = —2.5. 

Two-dimensional piecewise linear model 

The piecewise linear model is named for its piecewise linear nullclines in the phase space. The general form of the model we 
used in this paper is: 


4 

dt 

+ 

1 

II 

(16) 

dw 

dt 

= {k„v-w)/x„, 

(17) 


where /(v) is a piecewise linear function of v (membrane potential), w represents the recovery current, C is the membrane 
capacitance, = 5 is the time constant of w, and k„ = 0.45 is a coefficient. As a model demonstrate threshold phenomena in 
principle, the equations is dimensionless, so do the equations of 3D PWL model. 

We take /(v) as a piecewise linear function: 

{ kiv + bi, v<vi, 

kmV + bm, Vi<V<Vr, (18) 

krV + br, V > Vr, 

Slope rates of left, middle and right segments are represented as kj = —0.5, km = 0.5, and kr = —0.25, respectively, whereas 
the intercepts are bi = 0, bm = —1.5 and by = 17.25, and the parameters that separate the definition scales are v; = 1.5 and 
V, = 25.0. 
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According the derivation in Supplementary Note SI, the threshold voltages is: 


9 = — (w — i 


( 19 ) 


where kg is the slope rate of the threshold and is equal to 2k„Cl + y [k^x^ + C)^ — Ak„x„j is the and bg = 

{ie + bm){k„ — kg)/{k„ — km) is the intercept. The corresponding quasi separatrix in the plane of clamped voltage versus 
clamping time {Xc-Vc)'. 


9{Xc) 


-Ic_ 

Woe — bQ 

_ Tc 

(1—e ^■)k — kg 


( 20 ) 


Three-dimensional piecewise linear model 

We extended the two-dimensional piecewise linear model to three dimensional to show the separatrix in the three dimensional 
phase space: 


dv 

C— 

dt 

du 

dt 

dw 

dt 


f{v)-U-W + ie{t), 
{kuV u) / Xit^ 

{k„v-w)/x„, 


( 21 ) 

( 22 ) 

(23) 


where /(v) is a piecewise linear function of membrane potential v, /(v) has the same form as in formula ( fTSl l, u and w are the 
recovery variables, (t») is the time constant of u (w), and > 0 (k„ > 0) is a coefficient. 

Parameters of /(v) are the same as in the 2D PWL model, with the exception of k^ — 0.95, bm = —2.175, br = 57.825, and 
Sr = 50.0. The slope rates are ku = 0.45 and k^^. = 0.6. Time constants are x^ = 5.0 and x„ = 10.0. 

See the derivation in Supplementary Note S2, the explicit equation of the threshold can be written as: 


9 = —au — bw + 


1 -f akjt + bk^r 

ku T k\v kyu 


{bm T , 


(24) 


where a — 


and b = — 




with r,- and rj being the roots of Supplementary equation 


ri+rj+Xu(C-kmX„) ‘ ^ 

(S3). Although two threshold planes exist, only the plane shown in Supplementary Fig. S4 is located properly to serve as the 
threshold for this 3-dimensional PWL neuron. 


Classic Hodgkin-Huxley model 

The equations of the classic HH model ard^ 


tlv a 

Cm — = -gNani h{y - Eno) 

at 

-gKn\V-EK)-gL{V-EL) + ie{t). 
= {m^{V)-m)/Xm{V), 

^ = {h„{V)-h)/x,{V), 

^ = (n4v)-«)/T„(y), 


(25) 

(26) 

(27) 

(28) 


where T^(y) = I/ {ax{V) + j5x{V)) {x represents m, /z or n in the following equations and paragraphs) is the time constant of x, 
andx<x,(y) = ax(y)/ (o!j(y) -I- j3j:(y)) is the asymptotic value of x when time approaches -foo. The rate functions and are: 

dm =0.1 ) An = ^.Oe 18 ; 

«n=0.01jpiy^, ji„=O.U5e-lo- 

V 1 

Ctfi = 0.07£ ^ 1^1^ = go-V * 

eTC- + 1 
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Reversal potentials are Ef^a = 50mV, Ek = —llinV and Ei = —54AmV. 

The neuronal membrane potential is clamped to the set voltage by an ideal voltage-clamp protocol. Briefly, the voltage is set 
to be the clamped voltage while the other equations are allowed to evolve. By setting V to be constant in equations (|26ll-(|28ll, 
we obtain; 

x{Vc,t,xo) =Xoc{Vc) + {xo-Xoc{Vc))e . (29) 

This equation describes the asymptotic transient process of gating variable x when the neuron is clamped to the voltage Vc- 
Thus, we analytically calculated channel states xiVc, ic) for the initial states xq . x{Vc, Tc) is the new initial values to determine 
whether the neuron will firing or not following voltage clamping. Both the resting states and the subsequent evolution after 
voltage clamping were numerically calculated using the 4th-order Runge-Kutta method. The resting states were saved after 
200ms of running the model without any stimulation. After the voltage clamp, we did not provide other stimulus, thereby 
making the HH model an autonomous system again. 
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Figure SI . Threshold of step and ramp pulse determined hy separatrix in 2D-PWL model, (a) Trajectories (solid lines 
and dash-dotted lines) of step current injections with different amplitudes to the PWL model. Whether the change after a 
rectangular pulse injection can induce an AP is determined by whether the state have crossed the separatrix (red dash-dotted 
line). In this model, a minimum threshold amplitude (/,/„ ) must be exceed to have an transient ArHH (b) Trajectories of ramp 
current injection with different rate to PWL model. Whether AP will generate after the ramp current also determined by 
separatrix-crossing mechanism. The pluses’ duration of every solid line and the corresponding dash-dotted line with the same 
color differ 0.06. 
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Figure S2. An unstable separatrix of unstable node determines threshold voltages in 2D-PWL model, (a) Trajectories 
(solid lines with arrows) show the unstable separatrix of unstable node (the open circle at the bottom) does act as the threshold 
set of the system (the red dashed lines). Black dotted lines are the nullclines. (b) The local trajectories of the unstable node 
(expand v definition domain of middle segment in (a) to the whole v-axis). The parameters are: ki = —0.5, = 0.5, 

kr = —0.25 k„ = 0.6, vi = 1.5, Vr = 25.0, = 10.0, ie{t) = 0. 
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Figure S3. The time constant affect the distinctiveness of threshold voltage of PWL model. The state plane trajectories 
and time series of different initial voltages. The threshold phenomena become more distinct as increasing of the time constant 
of recovery variable w, which in subhgures are the following values: (a) & (b), = 2.0; (c) & (d), Th, = 5.0; (e) & (f) 

,Tw = 10.0 Other parameters are: ki = —0.5, k,n = 0.46, kr = —0.25, k^ = 0.35, v/ = 10, Vr = 80, ie{t) = 0. 
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Figure S4. Dynamic analysis of 3D-PWL model. Any pair equations of our 3D-PWL model define a nullcline in state 
space (black solid lines). Three nullclines intersect with each other in real or by elongating at three fixed points, however, only 
the one representing the resting potential is actually exists, the other two are virtual, i.e. not locate in their definition zone. The 
sepratrix separate the 3-dimensional space into two parts which defines the firing or non-firing region in the middle region. 
Arrows show the characteristic vectors, the red one is the one with negative characteristic eigenvalue. 
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Figure S5. Separatrix in 3D FHN model, (a) Fixed point (black point), nullcsurface (front surface), fold curve (black line) 
and the separatrix (surface behind the nullcsurface). (b). The canard trajectories passing through the points of fold curve 
determine the threshold. 
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Figure S6. Separatrix-crossing mechanism of HH model in m-h-n state space. The boundary of maximum voltages after 
voltage clamp mapped on voltage clamp surface in m-h-n state space numerically show the separatrix of HH model. Colors 
indicate different maximum voltages after voltage clamp. Yellow solid lines show voltage clamping processes at several 
voltages and the green dashed lines indicate same duration of voltage clamping, they are all marked with values at the end of 
the lines. On the separtrix, the instantaneous threshold (0) equal to clamping voltage 14. In the firing zone (hot colored region), 
we have Vc > 9 and in the blue region Vc < 9 when voltage clamping is off. 


23 , 


28 
















\ (ms) 


Figure S7. Points where clamped voltage and instantaneous threshold equal (i.e. Vc = 9 ) form threshold 
voltages/times hne. Left: The formation of minimum threshold voltage. Right: The formation of maximum clamping time 
which can induce AP. The criteria of AP is V,nax ^ —l5mV, so threshold points will form the —l5mV line in the contour map 
of maximum voltage (Figure S5). 
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Figure S8. Contour of maximum voltage after voltage clamp in Tc and 14 phase space. All contour lines consist of three 
segments: low voltage segment (indicating normal threshold), intermediate voltage segment and high voltage segment 
(overlapping with dV/dt = 0 at the moment the voltage clamping is off). 
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Figure S9. The deterministic version of threshold related to membrane potential and rising rate of depolarization. 

(a) Threshold for different clamped states. 0; {i = 0,1,2,3) are the corresponding threshold of i, {i = 0,1,2,3) respectively. 
The threshold increases with the membrane potential, (b) Threshold decreases with depolarizing rate of membrane potential 
induced by step current. The instantaneous current Io5{t) horizontally shifting of resting state have the lowest threshold for 
depolarizing current injection. From resting state point so is brought to threshold point 02 by step current with amplitude 0.164 
(Af = 11.67575), while the step current bring state from so to 0i having amplitude 0.2 (At = 4.7875). 
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Note SI: Derivation of the threshold curve of 2D-PWL model 

As shown in Fig. 3, v-nullcline have three segments and intersect with w-nullcline at three points (fixed points). This set of 
parameters makes the three fixed points, two stable spiral (focus) and a saddle, from left to right in the phase plane, respectively. 
The stable manifolds of saddle (separatrices) act as instantaneous threshold voltages in this PWL model. The separatrix s is a 
straight-line because the nullclines are linear. We assumed that the separatrix is w = kgr + fje, and we calculated the slope 

rate kg — 2k„C/ + C— ^ {km'^w + ~ and the w intercept bg = Then we obtained threshold 

voltages: 

0 = l(w-bg). (S30) 

kg 

In equation ( |S301 l, we replaced w with the transient value determined by the voltage clamping process: 

w = k„v+(wo — kt^v)e'^, (S31) 


then we obtained the corresponding quasi-separatrix in the clamped voltage versus clamping time (Zc-Vc) plane: 


9(Zc) = 


Woe — bg 

Tc * 

(1 —e^Tw)^ — 


The corresponding curve of dv/dt = 0 in the Zc — Vc plane also has a similar form. 


(S32) 
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Note S2: Derivation of the threshold plane of 3D-PWL model 

Shown in Supplementary Fig. S2, the piecewise linear function /(v) divides the whole space into three parts. For the above 
parameters, the left part has a real equilibrium, which is a stable node and represents resting state; the middle part has a virtual 
equilibrium in the left part, which is a saddle with its stable manifold consisting of thresholds; and the right part has a virtual 
equilibrium in the middle part, which is also a stable node and functions as a recovery mechanism. In order to show the 
dynamical properties well, we plot the nullclines (the intersect lines of nullcplanes dx/dt = 0(x = v, m,w)) only. 

The threshold set of our 3D-PWL model is a plane. See Figure 5 and S2, the plane is determined by stable (red arrow) 
and unstable (blue arrow) eigenvectors of a equilibrium (the blue circle). Each equation of 3D-PWL model defines a plane 
(three planes by the piecewise linear equation of dv/dt), whereas any satisfied pair of equations forms a nullcline in state space. 
These three nullclines intersect with each other at three fixed points; however, only the one that represents the resting potential 
is acmally exists, the other two are virtual, i.e., not located in their definition zone (see Figure S2). The saddle is 


(V/,M/,W/) = V/-(1,^„,^H,), 

where Vf = According to dynamical system theory, we obtain the following eigenvalues and 

characteristic matrix: 

A,- = 


'^u'^wC 


“U '’VV'-' 

_ / T tuC T {km kw) \ 

^ \ C'k T ^ k ^ ^ 

where r, (/ = 1,2,3) is the /-th root of the equation 

X* (C {Xii Tvv) km'^u'^w) ^ + (C T Tuik^ km) '^w{ku km)) 

+ {ku +k„- km)C^Tu'^l = 0 . 

Using the fixed-point coordinates and two eigenvectors, the threshold plane can be written 
kw'^uC (r/ f j km'^w)) ^/) k^X^ XwC{u uj) 

+ in in + nC) + T„ [rj + T„ (C - k„T„))) (w-Wf) =0, 

where i ^ j. So the explicit equation of the threshold can be written as 
0 = —au — bw + (1 +aku+bk„)vf, 
where a = , , ^ . and b = - ^ 

ri+rj+r„{C-k„T:„,) /:H.T„C(r,-+r;T„(C-i:„,T„.)) 


eigenvectors of the 

(533) 

(534) 


(535) 

using a point-norm form equation: 

(536) 

(537) 
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